In [7]:
from mpl_toolkits.basemap import Basemap
import pygrib as pb
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import glob
#Get data path
datapath = glob.glob("./21072000/M-*grb2")
print(datapath)
In [43]:
"""
利用basemap套件繪製圖形，在read_grib範例中，可以得知變數的位置，這邊一樣畫2米溫度跟10米風
需要的步驟，
1.讀取grib資料，取出所需資料，包含2米溫度、10米u風、10米v風、經緯度
2.投影相關設定值(lon_0, lat_0, lat_1, lat_2)
lon_0的key:LoVInDegrees
lat_0的key:LaDInDegrees
lat_1的key:Latin1InDegrees
lat_2的key:Latin2InDegrees
3.繪製圖形
"""
import sys
cnt = 0 # 只是為了抓取1測的固定值用
for grb in datapath:
if cnt == 0:
forget = pb.open(grb).select()[:][0]
analDate = forget.analDate.strftime("%y%m%d%H")
lon_0, lat_0 = forget["LoVInDegrees"], forget["LaDInDegrees"]
lat_1, lat_2 = forget["Latin1InDegrees"], forget["Latin2InDegrees"]
lats, lons = forget.latlons()
proj = Basemap(projection="lcc",resolution='l',lat_0 = lat_0 ,lon_0 = lon_0 ,lat_1 = lat_1, \
lat_2 = lat_2, llcrnrlat = lats[0,0], llcrnrlon = lons[0,0], urcrnrlat = lats[-1,-1], urcrnrlon = lons[-1,-1])
cnt = 1
basicvar = pb.open(grb).select()[:]
fcst = str(basicvar[0].forecastTime)
valDate = basicvar[0].validDate.strftime("%y%m%d%H")
t2 = basicvar[62]["values"] - 273.15
u10m_data = basicvar[66]["values"]
v10m_data = basicvar[67]["values"]
fig = plt.figure(figsize=(16,12))
proj.drawcoastlines(linewidth=0.9, color='slategrey')
ctf = proj.contourf(lons, lats, t2, cmap='coolwarm',latlon=True)
gap = 20 #沒有用gap的話，會很慘
proj.barbs(lons[::gap,::gap], lats[::gap,::gap], u10m_data[::gap,::gap], v10m_data[::gap,::gap] ,\
length=5,barbcolor="navy", latlon=True)
cbar = plt.colorbar(ctf,orientation='horizontal',pad=0.02)
cbar.set_label("$^\circ$C",size=14)
plt.title("Initial Time: {}".format(analDate), loc="left",fontsize=16)
plt.title("Forecast Time:{}".format(fcst),loc="center",fontsize=16)
plt.title("Valid Time:{}".format(valDate),loc="right",fontsize=16)
plt.suptitle("In-Fa Typhoon\n",y=0.93,fontsize=20)
plt.show()
plt.clf()
